Solution of the Percus-Yevick equation for hard discs 
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We solve the Percus-Yevick equation in two dimensions by reducing it to a set of simple integral 
equations. We numerically obtain both the pair correlation function and the equation of state for 
a hard disc fluid and find good agreement with available Monte-Carlo calculations. The present 
method of resolution may be generalized to any even dimension. 
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The description of hard-core fluids by approximate theories has attracted a lot of attention throughout the years 
The reason for this interest is twofold: firstly a solution of the full problem was, and still is, extremely difficult. 
Secondly, approximate methods have allowed for very good predictions in the low density phase. Among the most 
widely used approximations is the Percus-Yevick (PY) equation for d-dimensional hard spheres ■ 

Baxter developed a powerful method for solving the PY equation that works in odd dimensions [|[. Leutheusser 
used this method to reduce the nonlinear integral PY equation to a set of nonlinear algebraic equations of order d — 3 
(for d > 3) 0]. However, this set of equations can only be solved analytically for odd d < 7. In general, even the 
numerical solution of Leutheusser's equations is difficult to obtain for higher dimensions because a general solution 
to nonlinear algebraic equations of order larger than four is not available. Analytical results have been found for one 
i, three 0,0, five i and, most recently, seven dimensions [8|. 

In two dimensions an approximate numerical solution of the PY equation was found by Lado Q . Leutheusser (To| 
was able to fit many of Lado's results using an ansatz for the direct correlation function. Rosenfeld [ll[ generalized 
Leutheusser's Ansatz to higher dimensions and compared the results with the analytical results in three and five 
dimensions. All other results available in the literature are based on Molecular Dynamics (MD) or Monte Carlo (MC) 
methods 0, [H, Q [H M E3 ■ 

In this paper, we solve the PY equation for hard discs (d = 2). We develop a method that reduces the problem 
to a set of integral equations that are solved numerically without major difficulties. The originality of the present 
method is based on techniques borrowed from the resolution of crack problems [1 81 ] and uses some results from Baxter's 
| classical method [H, [j| . The main difference from the hard sphere case (and any odd dimension in general) is that the 
problem of finding the indirect correlation function and the direct correlation function are coupled. This means that 
the present analysis necessarily yields both correlation functions and therefore provides the equation of state. Note 
, that Lado [9( solved the problem using approximate integral equations, in contrast to our exact integral equations. 
The advantage of the current method is that one can improve the precision at will, and that it may be generalized to 
OO higher even dimensions. 
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. !-H , II. THE PY APPROXIMATION 

X: 

The pair correlation function g(r) is related to the direct correlation function c(r) through the Ornstein-Zernike 
equation by [l], Gi| 



h(r) = c(r)+p h(r')c(\r - r'\)dr' , (1) 
Jo 

where p is the particle number density and 

h(r) = g(r) - I , (2) 

is the indirect correlation function. The PY approximation is a closure relation of Eq. ([T]). For a hard-core pair 
interaction potential, this approximation reads 



I 



g(r) = h(r) + 1 - 0, r < 1 , (3) 
c(r) = 0, r > 1 . (4) 
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Here and elsewhere, we use the radius of the sphere as the unit of length. Thus in two dimensions, we have p = 4r)/ir 
where 77 is the packing fraction, and the space filling density corresponds to rj = 1. We define the two-dimensional 
Fourier transform 



F(q) =2tt / rJ (qr)F(r)dr , (5) 
Jo 

where Jo is the zeroth order Bessel function. The inverse Fourier transform is then given by 

POO 

F{r)^{2n)- 1 qJ Q (qr)F(q)dq . (6) 



In Fourier space Eq. ([I]) reads 

[1 - pc{q)\ [l + ph(q)] = 1 , (7) 
Finally, the static structure factor s(q) of wavenumber q is related to the pair correlation function through 

S (q) = l + ph(q) = — L_. (8) 

III. RESOLUTION 

Condition together with Eq. © imposes that c(q) can be written as [H| 

1 f 1 

c{q) = - cj)(t)sm.qtdt, (9) 



1 Jo 

where (f)(t) is a real function. Thus using Eq. ([B]), c(r) is expressed as [2C 



C M=r^M^. do) 

V ' J r v^^TT V ' 



The function c(r) does not diverge at r = 1 , thus the function cf>(t) satisfies 



Define 



4>{t) ~ 4 c(1 ^ as t -» I" . (11) 
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- A(q) = 1 - pc(g) = 1 - - / singidi . (12) 



s(q) ' ' "' ' '/ ./ , 

The function A(g) has the same properties as the corresponding function defined in the odd dimensional case [l[ : it 
has neither zeros nor poles on the real axis, since by definition s(q) has neither zeros nor poles for all g's, AW) = A(—q), 
and A(q) — > 1 as q — > 00. Therefore, in order to determine </>(£), we can use the formulation of Baxter J3|, which was 
developed to solve the PY equation in odd dimensions. We use the Wiener-Hopf method by defining [l|, y| 

A(q)=Q(q)Q(-q), (13) 

where Q(q) is an analytic function for 3(g) > 0. Following the same steps as in P, Q , it may be shown that Q(q) can 
be written as 

Q(q) = l-p [ Q(t)e tqt dt. (14) 
Jo 

Substituting (|13|14| into Eq. ^ gives 

- f 0(a) sin qsds = [ Q(s)e iqs ds+[ Q{s)e~ iqs ds ~ p ( ds [ ds'Q(s)Q(s')e iq{s - s " ) . (15) 
1 Jo Jo Jo Jo Jo 



Multiplying (fT5|) by exp(— iqt), with t > 0, and integrating with respect to q from — oo to oo gives 



fr(s)ds = 2Q(t) - 2p / Q(s)Q(s - t)ds . 
t Jt 

By setting t = 1 in ([IB]) we see that 

Q(1) = 0. 

Moreover, Eq. (|16p can be simplified further by differentiating with respect to t to give 

= -2Q'(t) + 2p [ Q'(s)Q(s - t)ds < t < 1 . 



(16) 



(17) 



(18) 



Therefore once the function Q(i) is determined, the functions 4>{t) and c(r) may be determined from Eq. (|18p and 
Eq. (fTO]) respectively. 

Now let us work on the function h(r). Since h(q) = h(—q), one can write without loss of generality 



qh(q) 



ip{t) sin qtdt , 



where tp(t) is a real function and h(r) is given in terms of ip(t) by 

ip(t) dt 



Substituting this form for h(r) in condition ([3|) yields 



(19) 



(20) 



< r < 1 



Vt 2 - r 2 2tt J x Vi 2 - r 2 2tt 

Eq. (|21[) is an Abel integral equation which can be inverted. It is easily shown that the inversion of the equation 

K ° <r(C) 



(21) 



yields 



d( 



d( = T((o) Co>0. 



c 



r(Co) dCo 



Therefore, using the change of variables ( = 1 — t 2 and ^ = 1 — r 2 , the inversion of Eq. (|2"Tj) gives 

-4* 



(22) 



(23) 



1 



Vs 2 - 1 , , . ds 



0<t<l 



(24) 



Eq. (|24|) is an integral equation fixing ?/>(£) for < t < 1 as function of i/;(f) for t > 1. 

Now, let us determine the relationship between Q(i) and VK*)- Putting together the results of Eqs. (|12I13I14I19[) in 
Eq. gives 



1 - P f Q{s)e tqs ds] [l + - / 
Jo 1 l 1 Jo 



1 + — / ip(s) sinqsds 
1 Jo 



Q(-q) 



(25) 



Multiplying Eq. ([23]) by exp(— ig<) with i > and integrating with respect to q from — oo to oo gives 0, [|| 



4Q(t)= dsip(s)(e(s + t) + e(s-t))- pi dsQ(s) ds'ip(s')(e(s' + 1 - s) + e(s' - t + s)) , (26) 
Jo Jo Jo 
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where e(x) — sign(a;). Using e(x) = — 1 + 20(a;), with Q(x) the Heaviside function, and differentiating Eq. ([26]) with 
respect to t gives 



2Q'(t)+i>(t)=p Q(s)[ip(t-s)-ip(s-t)]ds t>0. 



(27) 

Recalling that ip(t) is defined for t > and that Q(t) is defined for < t < 1, Eq. (|27|) can be split into two integral 
equations as follows 



2Q'(t) + V(i) = Q(s)i/)(t - a)ds - p J Q(s)ip(s - t)ds 

ip(t) = p [ Q{s)ip{t - s)ds for t > 1 . 
Jo 



for 



<t < 1, 



(28) 
(29) 



Our approach has reduced the PY problem for hard discs to the solution of the integral equations (|24I28I29[) with 
the additional boundary condition (TTTI) . We note that unlike the odd dimensional case Q, here one cannot separate 
the problem of finding the direct correlation function c(r) from that of finding the pair correlation function g(r). This 
is because the behavior of ip(t) for < t < 1 is related to the behavior of i/)(t) for t > 1 through Eq. (|24p. Although 
we were unable to find an analytical solution valid for all p, the numerical solution of these equations can be easily 
implemented. Before dealing with the numerical analysis, let us first consider the equation of state in the present 
formulation of the problem. 

There are two methods used to calculate the equation of state when the radial distribution function, g(r), is known. 
Without the assumptions made in deriving the PY equation, these two methods would yield the same equation of 
state. The difference in the equations of state calculated using these two methods therefore provides an estimation of 
the error made by using the PY approximation. The first equation of state is derived from the virial theorem and is 
given by [l[ 



(3P 



P+\ P 2 g(l + ). 



Using g(l + ) = — c(l ) 1] and Eqs. (|llll8p . the "virial" equation of state becomes 

PP (v) =p+^-p 2 lim y/1 - t 2 Q'(t) . 
4 t->i- 

The second method uses the isothermal compressibility kt which is given by [l[ 

dp 



pp- 1 ^ = 13- 1 



ap(c) 



= s(q - 0) 



(30) 



(31) 



(32) 



Using Eq. (fT2|) and replacing <p(t) with Q(t) by using Eq. (fl~8|) . one has 



1 



l + 2p 



p / Q{s)Q(s-t)ds-Q(t) 



dt. 



Integrating (f3"2"]) and substituting for s(0) 1 from (f3"3"]) we find the "compressibility" equation of state 



/3P 



rp ( rl 
p + 2 2p' 

Jo (JO 



// / Q(s)Q(s-t)ds-Q(t) 



dt^dp' . 



(33) 



(34) 



IV. NUMERICAL PROCEDURE 



We found it simplest to implement an iterative procedure for solving the set of integral equations (|24l28l29p . First, 
note that for p = 0, the exact solution is given by 



m 

Q(t) 



-At 



6(1 -t) 



-2y/l - i 2 6(l - t) 



(35) 
(36) 
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and that for any p, one can write the functions Q(t) and ip(t) as power series in p 

m = = ^ ~ J ^p^(')(t) + ^pV(')ft)e(t-l)8(» + l-t). (37) 



Q(t) = -2 vT^t* ( 38 ) 



i=0 

with the definition 



= g(°)(t) =1, (39) 

V> (0) ft) = 0. (40) 



Using this power series representation, Eq. (|24|) yields 



= i>l, (41) 



and Eq. (|28|) becomes 



(l-t 2 )g'( J + 1 )(i)-<< Z ( l + 1 )(t)+t*( l+1 )(i) = -2VT T ^ / ^ ' S)Vl Z y]gW(a)$^-fc)(|t-a|)ffa l >0. (42) 

Jo y/l-(t-s) 2 ^ 

The boundary condition q^(l) = for Vi > is necessary to avoid a singularity in q^ l \t) as t — > 1. Finally, 

Eq. {22]) gives 



t-i VI - s 2 7i 



1 < 



fc=0 

ij (l+1) (t) = -2^2 [ y/l-(t- s) 2 q {k) (t - s)^ k) {s)ds, t>2, (44) 



In devising a successful numerical scheme for this problem, it is important to note that several of the integrands 
have mild singularities, for example at s — 1 in (|4"3"| and ats = l,t=lin (14"T]) . The singularity in (I43[) may be dealt 
with simply by subtraction. To deal with the singularity in (|41[) , we note that we may integrate (|41| once by parts 
and write 



ruu, n o f' l+1 /s 2 - 1 d fip^(s)\ ds f l+1 pr, -df^(s)\ds 



V s J 2tt 



(45) 



This obviates the need to subtract the singularity and also suggests that we may write 

*(*) ft) = + $(0 Vl-t 2 + g {V > ft), (46) 

with ,9 (l) (l) = 0, and 

"■af-^(*)i-^ 

Examining the equation for q^ l+1 \t), (|4"2")l . we see that the form assumed for (|4"6"]l suggests that 

gW ft) = a « + 7 W Vl-< 2 + f {l) ft) . (49) 
Here, the function /Wft) satisfies / w (l) = and 



(1 - f 2 )/'( 4+1 >ft) + i[.g (l+1) ft) - / (4+1) ft)] = 2-n/I -i 2 [iA (4+1) (l) - A^ft)] , (50) 
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where 

A^\t)^±f {t ~ S) f^Z ^\ S )^- k H\t-s\) ds (51) 

The constant may be determined by considering the limit t — > 1, which gives 

7 W = !<5« + (52) 

We solve the system of equations by discretizing in space using steps of size A. Integrals are evaluated using the 
trapezoidal rule and derivatives are calculated using forward-differencing, both of which are first-order accurate in 
space. Iterations proceed from i = using the values \E , ^°- ) and q(°\ given in (|40|) . to determine ipW(t). Using ^^(t) 
in the discretized versions of (|47|) and (j48| allows us to determine ^^(t). Substituting these into the discretized 
(|50[) . q^\t) may be determined. The whole process is then iterated N times, corresponding to determining the first 
N terms in the series for the variables q(t) and ijj(t). The results presented here typically have N = 20. 



V. NUMERICAL RESULTS 



Let us begin by presenting our results for the virial coefficients and the equation of state. The virial coefficients, 
Bi, are defined by 

00 

(3P = J2BiP i . (53) 

It is well-known that B\ — 1 and B2 = 7r/2 are exact and are recovered within the PY approximation. The higher 
order virial coefficients can be calculated from the solution of the PY equation using the expressions in (|3"Tj) and (I34|) . 
For example, using Eqs. (|38I31[) yields 

flW = ^>2, (54) 

(c) 

Thanks to the iterative procedure presented above, these coefficients as well as the B\ , which are derived from the 
isothermal compressibility (|34p , are directly given by the numerical resolution of the problem. 

Our numerical solution of the PY equation and calculation of the virial coefficients are only first order accurate in 
space and so we expect to accumulate errors at 0(A). To counter this, we calculated each of the virial coefficients 
for several values of A and then extrapolated linearly its value at A = 0. For each virial coefficient the correlation 
coefficient was calculated to be < —0.999 confirming that the expected linear dependence of Bi on A is indeed 
observed. The computed values of the first twenty coefficients for the series corresponding to f3P v and (3P C are 
presented in Table HI The values of the first ten virial coefficients for hard discs were determined recently using a "hit 
or miss" Monte-Carlo integration algorithm in [l7| and are reproduced in Table U for comparison. 

As is clear from the definition (|53p . the virial coefficients are important for determining the equation of state. The 
value of the pressure, (3P, is therefore of considerable interest. Fig. [T] shows the dependence of j3P on the packing 
fraction of the disks, 77. This shows that there is good agreement between our results and previous ones using a direct 
approximate numerical resolution of PY equation |9|. Furthermore, this agreement is within the 2% error estimate 
given in 0. 

It is clear from Fig. [T] that there is a divergence in (3P for relatively large 77. The question of interest here is: at 
what density does this divergence occur? We are thus interested in determining the radius of convergence of the series 
for pP. To this end, we plot in Fig. [5] the ratio between successive virial coefficients given in Table fl] In particular, 
we observe that as i — > 00 

B i+ 1 7T 

Using d'Alembert's ratio test, this observation suggests that the series for j3P will converge absolutely provided that 
p < 4/7r or, equivalently, r/ < ?/ c = 1, the space filling density. This result shows that, similarly to the case of hard 
spheres (d = 3), the PY equation for hard discs predicts no phase transition at intermediate density and thus fails at 
high densities. 

Finally, the correlation function g(r) may be calculated numerically from the values of tp^ . Some typical results are 
shown in Fig. [3] and compared with the available Monte Carlo results [H, [l5[ ■ As can be seen, the overall agreement 
is very good, except in the vicinity of r = 1. Also note that our solution of the PY equation in Fig.[3[a) agrees better 
with the Monte Carlo results than a previous numerical solution of the PY equation presented in [13| . 
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B (MC) 


3 


1.930 


1.930 


1.930 


4 


1.941 


2.084 


2.063 


5 


1.795 


2.090 


2.031 


6 


1.594 


1.999 


1.902 


7 


1.381 


1.852 


1.726 


8 


1.177 


1.677 


1.534 


9 


0.992 


1.492 


1.342 


10 


0.829 


1.309 


1.162 


11 


0.688 


1.136 


— 


12 


0.567 


0.977 


— 


13 


0.466 


0.833 


— 


14 


0.381 


0.707 


— 


15 


0.311 


0.596 


— 


16 


0.253 


0.500 


— 


17 


0.205 


0.418 




18 


0.16(6) 


0.348 




19 


0.13(4) 


0.28(9) 




20 


0.10(7) 


0.23(9) 





TABLE I: Numerical values of the first twenty virial coefficients. The B\ are the results from Monte Carlo calculations 
presented in [1?]]. Bj and b\°' are the values found from the solution to the PY equation using Eqs. (|54[1 and (|34l) respectively. 

25 
20 

r— I 

^10 

5 



"0 0.2 0.4 0.6 0.8 

FIG. 1: Reduced pressure^P/p — 1 as a function of the disk packing fraction, rj. The results of molecular dynamics simulations 
reported previously [13 . Il4 ] are shown as + along with the results obtained from an earlier solution of the PY equation 
shown as x . The two solid curves show the results of our calculations using the first 20 virial coefficients, with the labels 
corresponding to the equation of state used to calculate them. 




VI. DISCUSSION 



In this paper we developed a semi-analytic method to solve the PY equation for hard discs. At the heart of this 
approach is a reduction of the PY equation to a set of integral equations for two auxiliary functions Q{s) and ip(s) 
as given by Eqs. (|24l28l29p . The correlation functions and the equation of state can be determined easily from these 
auxiliary functions. We suggest an efficient iterative numerical method to solve these equations and determine the 
auxiliary functions. Using this method we are able to determine the values of the virial coefficients within the PY 
approximation. Furthermore, it allows for a comparison with the available 10 virial terms of the full problem p7| . 
We also obtain results for the pair correlation function which compare well with the available MC calculations. 

An advantage of this method is that it allows, in principle, calculations to arbitrary precision, and it could be 
interesting to obtain more virial coefficients by so doing. It could also be interesting to generalize this approach to 
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FIG. 2: The ratio between successive virial coefficients Bi/Bi+i as a function of i. The ratios determined from the solution to 
the PY equation are denoted by + ((3P V ) and x (/3P C ). These results sandwich the virial coefficients determined from Monte 
Carlo calculations given in [13], which are represented by O- The dashed line represents B t = 4B i+ i/n, which appears to be 
an asymptote as i — » oo. 




FIG. 3: The correlation function g(r) computed from the PY equation (curves) and from different Monte Carlo simulations 
(crosses), (a) p = 0.462: Our solution of the PY equation (solid curve) compared to the MC results of Chae et al. [H| (crosses) 
and their solution of the PY equation (dashed curve), (b) p = 0.794: Our solution of the PY equation (solid curve) compared 
to the MC results of Wood fl5l ] (crosses). Solid curves were calculated using the first 50 terms of the series with A = 0.00125. 



polydisperse mixtures [2l| and to higher even dimensions. 
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